library(tidyverse)
library(plotly)
library(sf)
library(tigris)
library(leaflet)
library(censusapi)
library(mapview)
library(png)
library(raster)

Sys.setenv(CENSUS_KEY="c8aa67e4086b4b5ce3a8717f59faa9a28f611dab")

The data comes from http://www.jnoble.org/ECOSTRESS/2020-08-03/baseMap1/index.html

ECOSTRESS Land Surface Temperature (LST) Silicon Valley, CA, USA – 2020.08.03, 9:42 pm PT Data: Glynn Hulley, ECOSTRESS team, NASA JPL

Load thermal data

From layers.js in the .zip file of “Land parcels and buildings, 1:25,000 (OL web map, 230 MB), (.zip, 44 MB)” in http://www.jnoble.org/ECOSTRESS/index.html, I noticed the following details about the PNG:

projection: ‘EPSG:3857’,

imageExtent: [-13616154.835870, 4484784.913060, -13582758.988632, 4509845.252541]

So I looked up some techniques online, played around a bit, and was able to load in the PNG as an array, line it up with pixel coordinates in the right CRS, and then create a raster. Note I’m only using ecostress[ , , 2] which is only the “G” of the “RGB” data, which looks like it’s stored as fractions (we’ll need to figure out how to convert this back to thermal information, which should be something we can figure out from the literature on ecostress). There’s a flip() I need to do to get it in the right orientation – not exactly sure why but it looks like it’s in the right orientation at the end of this process. From here we can then do some work with the raster package, similar to how we intersect raster images with building footprints in this flood analysis: https://github.com/stanfordfuturebay/surge20/blob/master/SURF.Rmd#L80.

ecostress <- readPNG("G:/Shared drives/SFBI/CRNFO/ECOSTRESS_2020-08-03_SV_webmap1/ECOSTRESS_2020-08-03_SV_webmap1/layers/ECOSTRESSLSTC_2.png")

dat1 <- NULL

dat1$x <- seq(
  from = -13616154.835870,
  to = -13582758.988632,
  length.out = 932)

dat1$y <- seq(
  from = 4484784.913060,
  to = 4509845.252541,
  length.out = 700)

dat1$z <- ecostress[ , , 2] %>% t()

final <- raster(dat1) %>% 
  flip(direction = 'y')

crs(final) <- CRS('+init=EPSG:3857')

mapview(final)